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There has been recent progress on inferring the structure of interactions in complex networks 
when they are in stationary states satisfying detailed balance, but little has been done for non- 
equilibrium systems. Here we introduce an approach to this problem, considering, as an example, 
the question of recovering the interactions in an asymmetrically-coupled, synchronously-updated 
Sherrington-Kirkpatrick model. We derive an exact iterative inversion algorithm and develop effi- 
cient approximations based on dynamical mean-field and Thouless- Anderson-Palmer equations that 
express the interactions in terms of equal-time and one time step-delayed correlation functions. 
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Introduction. — Finding the connectivity in complex 
networks is crucial for understanding how they oper- 
ate. Gene and multi-electrode microarrays have recently 
made the type of data required for this purpose available. 
What is needed now is appropriate theoretical tools for 
analyzing these data and extracting the connectivity. 

In much recent work on this subject [H-Qi the prob- 
lem has been posed as that of inferring the parameters 
of a stationary Gibbs distribution modeling the system. 
While satisfied in many applications, the assumption of 
Gibbs equilibrium is unlikely to hold for many biological 
systems since they are usually driven by time-dependent 
external fields, their interactions may not satisfy detailed 
balance, or they may only be observed while the tran- 
sients dominate the dynamics. Applying the equilibrium 
approach to such cases usually yields effective interac- 
tions that do not bear an obvious relationship to the real 
ones Kinetic and nonequilibrium models provide a 
much richer platform for studying such systems 0-01 ■ 

Whereas for equilibrium models the development of 
systematic mean field inference methods Q has led to 
great practical and conceptual advancements, a mean 
field theory for nonequilibrium network reconstruction 
is still lacking. In this paper, we show how a mean field 
theory for inference can also be developed for a nonequi- 
librium system. We consider this problem for a particular 
simple nonequilibrium model: a kinetic Ising model with 
random asymmetric interactions (Jji independent of Jij), 
in an external field which may be time-dependent. This is 
a discrete-time, synchronously updated model composed 
of TV spins Si — ±1 with transition probability 



Pr(s(t+l)|s(t)) = n 



2cosh(6'j(t)) 



(1) 



where 9i{t) = hiit) +^jJijSj{t). The couplings Jij 
are independent Gaussian variables with variance /N . 
This model can be readily applied to time-binned neu- 
ral data, where t labels the bins, and Si{t) = ±1 repre- 
sents a spike or no spike by neuron i in bin t The 
temperature has been set equal to 1, since any effects of 



changing the temperature can be realized by changing the 
coupling parameter g and the field strengths. Even for 
time-independent field and in a steady state, this system 
is not in a Gibbs equilibrium [lo| . However, we show 
that, like its equilibrium counterpart, the nonequilib- 
rium inverse problem for this model can be solved using 
a gradient descent method and also via systematic ap- 
proximate inferences derived using dynamical versions of 
naive mean-field (nMF) and Thouless- Anderson-Palmer 
(TAP) equations. We show that for both the stationary 
and nonstationary systems these methods provide effi- 
cient reconstruction of interactions. We also analytically 
quantify their errors. 

Exact, nMF and TAP learning rules. — Sup- 
pose that we have observed R realizations of duration 
L time steps of the process in ([T]). We denote the ob- 
served state of the system at time t of realization r by 
s^[t) = {s\{t), • • • , s]y(t)}. To find the couplings and ex- 
ternal fields, we maximize the likelihood of the observed 
states under the model ([l}. This maximization can be 
done using an iterative algorithm, analogous to Boltz- 
mann learning for the equilibrium model: starting from 
an initial set of couplings and fields, one adjusts them it- 
eratively by steps of sizes 5hi = rih§^ and SJij = VJ'^j^j 
L being the log-likelihood. The learning steps thus are 

^h,{i) - i-^h{{s,{t + 1)),. - (tanh[0,(i))]),]} (2a) 
5J^, = ??j{(s,:(i + l)sj(t)) - (tanh[0,(t)]sj(t))} (2b) 

where 77/1 and r\j are learning rates. Here and in what 
follows (• • • )ri (■ • • ) represent averaging over repeats, and 
both repeats and time, respectively. An overline, instead, 
will indicate averaging over the spins. One can think of 
Eq. (j2bp as performing a logistic regression to explain 
one-step separated correlations. This is similar to what 
is proposed in [13[ as an approximation for inferring the 
connectivity in an equilibrium Ising model. 

Since performing the steps in this algorithm does not 
require Monte Carlo runs, it is faster than the equilibrium 
Boltzmann learning. However, two factors still make this 
algorithm slow for large systems and/or data sets, war- 
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ranting the development of fast approximations. First, 
([2]) is still an iterative algorithm which could take a long 
time to converge if not provided with a good initial condi- 
tion and learning rates. Second, at each step the averages 
on the right hand side of ^ should be calculated from 
the data de novo, given the adjusted parameters. 

Two fast approximations, nMF and TAP learning 
rules, are derived and studied below. To implement them 
in the stationary case, one first uses the data to calcu- 
late the one-step delayed and equal time correlations, 
Dij = {Ssi{t + l)6sj{t)) and Cij = {Ssi{t)Ssj{t)) , where 
nT-i = (si) and Ssi ~ Si — rrii. The approximations are 

jnMF/TAP _ y^nMF/TAP^lQ(--l 

where Aff ^ = (1 -m?)^^, Ajf^^ = A'lJ^^{l-F,) and E, 
is the root of the cubic equation ([6]) below. In the non- 
stationary case too, similar learning rules can be derived 
as shown later in the paper. 

Derivation of nMF and TAP inversion. — For 
simplicity, we consider first the stationary case, for which 
the sequence index r is superfluous, as averaging over 
time and repeats would be equivalent. We start with 
the maximum likelihood conditions, i.e. Shi = SJij = 
in ([2]). Using the nMF equations rrii = tanh(/ii + 
writing the Si in ^ as nii + Ssj, we 
expand the tanh in the Ssi. The first nonzero term gives 

{Ss,it + l)Ss,it)) il-mf)Y,J-nSskit)Ss,it)). (4) 

k 

which can be written as ([3]) for the nMF case. 

To get the TAP inversion formula, we start instead 
by assuming that the mi satisfy the TAP equations 
= tanh[/i, + Jlff^'mi. - " ^1)], 

which take into account the Onsager reaction term. Kap- 
pen and Spanjers Q proved that the TAP equations, al- 
though usually derived for the equilibrium (symmetric- 
J) SK model, also hold for the asynchronously updated, 
asymmetric-J model in a stationary state. We have veri- 
fied that they are also valid in our synchronously-updated 
model 0. We again write Si = + Ssi, expand- 
ing the tanh to third order in powers of J^^Ssk + 
™j ~ "^fc)- Keeping terms up to order p"^ 

leads to D = ATAPjtap^ ^ ^1^^,^.^ 

A^AP = ^=^^^[1 - (1 - m?) ^( JTAP)|(1 _ mf)]. (5) 

I 

These equations cannot be solved directly as in the nMF 
case because A"^"^^ depends on J'^^^. However, one 
can derive a cubic equation for the quantities Fi — 
(l-mf)E,(J^^")i(l-mf): 

- FO^ = (1 - m?) ^(J"^^)?,(l - m'^). (6) 
j 

This determines Ajf"^ = ^"'^^(l - F,), yielding © for 
the TAP case. The relevant root of ([H]) is the smallest one 



(the one approaching zero as g ^ 0). This root cannot 
exceed 1/3, restricting this technique to weak couplings. 

For both nMF and TAP reconstruction, the external 
fields hi can also be found by solving the respective mag- 
netization equations after the Jy have been obtained, 
just as in the equilibrium problem 0- 

Performance of the algorithms. — We have veri- 
fied that the algorithm ([2]) recovers the couplings of an 
asymmetric SK model exactly in the limit of L — !■ oo, for 
a wide range of coupling strengths g, external fields and 
system sizes. The mean square error, Cexact, is in gen- 
eral proportional to 1/L, and in the weak-coupling limit 
a quadratic expansion of log-likelihood yields 

Ccxact = Sjfj EE {Jij - = _ ^2^2.' 

where JijiJi^'^) are the inferred (true) couplings. 

We find that the nMF algorithm leads to an error, cmf, 
of the form ecxact + £5^^IF' where is independent of 

L and proportional to Thus, for data sets of length 

L <^ L* ~ I/cJJmf °c N, nMF does almost as well as the 
exact algorithm. Furthermore, the larger the network, 
the better nMF does. The errors for the exact and nMF 
algorithms vs L are shown in Fig. 
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FIG. 1. (Color online) Performance of the algorithms. Exact 
and nMF (a), and the TAP (b) erros are shown vs data length 
L for g = 0.1 (blue stars), 0.12 (magenta crosses), 0.14 (red 
circles) and 0.16 (black x), all for TV = 20 and zero external 
field. Theoretical predictions are the solid lines. 

For weak coupling, we can calculate the asymptotic 
nMF error, e^mF' sualytically as follows. We present the 
zero-field case here for simplicity. We expand the tanh 
in the max-likelihood equation to third order, giving 

E^in — ^ ^ Jiki^k^n) 3 ^ ^ Jik Jil J hn {^k^l^m^n) ~^ ' ' ' • 
k klm 

(8) 

Correlations here are at equal times, except for Din- The 
dominant contributions in the sum over fc, /, m arc those 
with k = I, I = m and m = k. Multiplying on the right by 
(C^^)nj, summing over n and using ([3]) for nMF, yields 

Jij — Jij — ^ ^ JikJij, (9) 

k 

with corrections of relative order Eq. © also yields 
the TAP-approximation couplings found above, showing 
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that the TAP reconstruction indeed corrects the leading 
MF errors. To leading order the sum on k is just g^, and 
the asymptotic mean square MF error is 



CMF = (^..-^f^)^-^- (10) 

The solid curves in Fig. [TJi are 1/L + g^/N; the fit is 
evidently good. As shown in Fig. SI nMF exhibits 
a systematic error by underestimating the magnitude of 
the couplings. The factor 1 — i^i in TAP formula cor- 
rects for this to relative order g'^. Thus, when one is in- 
terested only in the presence or absence of connections, 
there would be little difference between nMF and TAP. 

The error for the TAP reconstruction is much lower 
than that of the nMF one and reaches its minimum at 
much larger L: for = 20 and the coupling strengths 
we studied, we had to go to L ~ 10^ to see the error flat- 
ten (Fig.[T}D). To calculate the asymptotic reconstruction 
error for TAP, we expand the tanh to 5th order and pro- 
ceed to evaluate the averages as we did for nMF. The 
nMF error terms analyzed above are compensated for by 
the TAP equations, as A^ — >■ oo, leading to an asymp- 
totic e^p = 4g^°/N. For A^ > 1/g^ this is the leading 
term in the asymptotic TAP error. Outside this regime, 
a finite-size effect should also be taken into account. This 
is because in making that TAP correction, the term in (jSj 
with k ~ I = m has been counted three times in obtaining 
© instead of once. The mean square error that results 
from this overcounting is (2/3)2 = {2Qg^) / {?,N^) and 
should be added to the Ag^'^ /N term. 

Non-stationary case. — The magnetizations, 
niiit) = {sl{t))r-, are now time-dependent and, for nMF, 
solve 

m,{t 1) = tanh[/i,(t) + ^ J^J^^mj{t)]. (11) 

j 

We have also proved Q that the TAP equations hold 
even in a nonstationary state, in the form 

m,{t + 1) = tanh[/i,(i) + ^ J^^^mj{t) 

j 

- m,{t + 1) JT^P)5(1 - (12) 

j 

Thus, we can extend both our inversion algorithms to 
nonstationary systems, as we show in the following. 

We start by defining time-dependent correlation ma- 
trices Dij{t) = {5sl{t + l)5sj(t))r and C,j{t) = 
{Ssl{t)Ssj{t))r- For nMF, using the same procedure that 
lead to we find 

(A, - E J^''^^^ - ^"'(^ + l))Cfc. W)*- (13) 

k 

One can still solve for J by simple matrix algebra: 

jnMF ^j2{D,kimiB'^'^r%j, (14) 

k 



where B^^^ = ((1 - mf{t + l))Ckj{t))f The problem is 
more complex than the stationary one only because one 
has to invert a different matrix B*^'^ for each i. 

For TAP, analogously to the stationary case, the B*^*) 
acquire an extra factor inside the time average: 

Bis = - ™'(^ + - m))CkAt))u (15a) 
F,it) = (1 - mUt + 1)) - "^'(^))- (15b) 

Exact TAP inversion requires iterative solution of (|14p , 
with J^^P instead of J-j'^^, together with We have 
found, however, that effective reconstruction is still pos- 
sible under the simplifying approximation that Fi{t) in 
Eq. (|15ap can be represented by its temporal mean. In 
this case, F,; = {Fi{t))t solves the cubic equation 

F,(l - E, f = E(^"''^)| «l - "^'(^ + !))(! - "^?(*)))*- 

3 

Solving it and using it in Eq. (jl5ap . one can calculate 
JTAP = J^J.MF/(1 - F,). Similar to the stationary case, 
after inferring the couplings, one can use the forward dy- 
namical nMF and TAP equations Eqns. (dH and 
to infer the time-varying external field. The result of 
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FIG. 2. (Color online) Inference in the nonstationary case, 
(a) Couplings of a network of A'^ = 20 driven by a sinusoidal 
external field inferred using the nonstationary nMF, and (b) 
the stationary nMF. (c) Two periods of the external field (thin 
blue full curve) and its reconstruction using the nonstationary 
nMF couplings (red dashed curve) and stationary nMF (thick 
black full curve). 

reconstructing the couplings of a network driven by a 
common sinusoidal external field to all spins is shown in 
Fig. [5] Fig. [2^ shows how well the couplings are inferred 
by nonstationary MF using L = 10^ and R = 100. Non- 
stationary TAP couplings (not shown) have a lower mean 
squared error: 6.7 x 10~^ versus 10~^ for nMF. In Fig. 
[2t), we also plot the couplings inferred using stationary 
nMF inversion for each of the 100 repeats and averaging 
over them. Not surprisingly, the stationary nMF per- 
forms poorly on this nonstationary data. Importantly, 
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there is a systematic overestimation of the coupUngs in 
this case, because the stationary method accounts for cor- 
relations induced by the common, time-varying external 
field through adjusting the couplings. Correspondingly, 
if one uses the couplings inferred by stationary nMF in 
(fTTjl to infer hi{t), the amplitude of this field is underes- 
timated, while the use of nonstationary nMF couplings 
yields a very good reconstruction of hi{t); see Fig. [2t. 

Discussion. — We have shown how to infer interac- 
tions in a simple but nontrivial nonequilibrium system: 
a kinetic Ising model with random and potentially asym- 
metric interactions. The model is the maximum entropy 
model for each time step, given mean magnetizations and 
one step separated correlations. We have described both 
an exact iterative algorithm and two approximate ones, 
based on dynamical nMF and TAP equations, which are 
correct up to corrections of order 1/iV. We calculated 
analytically the errors of these approximations for weak 
coupling. The method shows particular promise when 
applied to nonstationary states, where it separates true 
interactions from the apparent ones found by applying a 
stationary theory to a nonstationary state. 



a) 

8 



X 10 


b) 

6 


X 10 














4 








II-- 


2 






-0.2 -0.1 0.1 0.2 ^ 

.TAP 


-0-2 -0-\tap<: 


0.1 0.2 



FIG. 3. (Color online) Finding connections in a cortical net- 
work model, (a) The histogram of the couplings inferred us- 
ing the stationary nonequilibrium TAP for pairs of neurons 
that were connected (blue full bars), and those that were not 
(red empty bars). The separation between the histograms 
shows that one can use the TAP approximation to separate 
connected and disconnected pairs, (b) same as (a) for equi- 
librium TAP. 

A kinetic Ising model will show an intrinsic error when 
applied to data from a different kind of system. However, 



resulting distributions are completely overlapping when 
and equilibrium TAP learning is used. When using a 
model like (H]) to infer connectivity in a system with a 
different dynamics, or when faced with data limitation, 
including prior knowledge about the network could be 
very beneficial. In particular, taking into account spar- 
sity of the connections via a l-l regularizer added to the 
likelihood has been shown to be very useful [l^. It is 
easy to show that adding an l-l regularizer to the likeli- 
hood of the data under ([T|) would modify ^ by adding a 
term proportional to A"'^^/"'"^^ sgn(J)C^^ to the right 
hand side. How this improves inferring connections in 
biological networks will be discussed elsewhere. 

A simple extension of ([TJ is its continuous time version. 
As shown in 15|, for this model, too, a mean field theory 
can be developed using the approach presented here. In 
other recent kinetic approaches to problems like this, the 
equilibrium maximum-entropy approach [l[ is extended 
to include non-equal-time correlations [Hj and an approx- 
imate scheme for fittingan integrate-and-fire network to 
data was developed in ^ . There has also been work , 
closely connected to ([I}, in which Si(i-I-I) depends on lin- 
ear combinations of h{t') and s(t'), for t' < t. Given the 
advantage of these nonequilibrium models over the equi- 
librium ones for describing spike train statistics, a mean- 
field theory for inferring their parameters would be of 
great theoretical and practical benefit. For such models, 
we expect that it will be possible to use the techqniues in 
[i or m to derive dynamical nMF and TAP equa- 
tions. Employing the approach developed here one can 
then build approximate mean field inversion techniques 
based on these equations. 
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